library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(lme4)

options(scipen=999)
df = read.csv('avg_tensor_by_roi_wide.csv',
              colClasses = c('subject' = 'factor',
                             'site' = 'factor',
                             'visit' = 'factor')) %>% 
  rename(scan = visit)
outcomes = df %>% 
  select(where(is.numeric)) %>% 
  colnames

# replace ATROPOS measures with NA for select images (segmentation failed)

atropos_cols <- df %>% 
  select(contains('ATROPOS')) %>% 
  colnames()
df[(df$subject == '04001' & df$site == 'NIH' & df$scan == '01') |
   (df$subject == '01003' & df$site == 'NIH' & df$scan == '01') |
   (df$subject == '03002' & df$site == 'NIH' & df$scan == '02') |
   (df$subject == '04003' & df$site == 'NIH' & df$scan == '01') |
   (df$subject == '04001' & df$site == 'BWH' & df$scan == '02') |
   (df$subject == '03001' & df$site == 'BWH' & df$scan == '01'), atropos_cols] <- NA

Summary

The purpose of the DTIMSCAMRAS study is to test for site effects on DTI metrics of different brain regions in the context of a traveling subjects design with a harmonized imaging protocol. Participants (n=11) with stable Multiple Sclerosis (MS) were scanned in 4 different locations: National Institute of Health (NIH), Brigham and Women’s Hospital (BWH), Johns Hopkins University (Hopkins), and The University of Pennsylvania (Penn). Two scans were collected per site visit, and participants were re-positioned between scans. The modalities that were collected include T1s (MPRAGE), FLAIR, and Diffusion-weighted images (DWI).

Preprocessing

To preprocess imaging data, qsiprep was used for bias correction, skull-stripping and co-registration of the T1s and DWIs, and specialized denoising and artifact correction for the DWIs. Then, DTIFIT was used to compute the following scalar maps: * Fractional Anisotropy * Mean Diffusivity * Axial Diffusivity * Radial Diffusivity

Concurrently, several segmentation pipelines were used on the T1s (+ FLAIRs in the case of MIMOSA) with the purpose of defining particular regions of interest (ROIs):

  • White Matter and Gray Matter Segmentation methods:
    • ANTs ATROPOS
    • FSL FAST
    • Multi-Atlas Segmentation with Joint Label Fusion using the MUSE templates (JLF WM and JLF GM)
  • Thalamus Segmentation
    • FSL FIRST
    • Multi-Atlas Segmentation with Joint Label Fusion using the OASIS templates (JLF THALAMUS)
  • Lesion Segmentation methods:
    • MIMOSA

The segmentation results were used to average the scalar maps across all voxels belonging to each of the ROIs. The 36 resulting metrics make up the set of outcome variables which are the focus of the site effects test. Below, the outcome variables are listed.

data.frame(OUTCOME = outcomes) %>% 
  separate(OUTCOME, c("SEGMENTATION", "ROI", "SCALAR"), remove = FALSE)

Note the general format; for example, the ATROPOS_GM_AD variable denotes the average Axial Diffusivity (AD) across the Gray Matter (GM) regions defined by the Atropos Segmentation.

Testing site effects

A permutation-based F-test was used to test for site effects in the absence of parametric assumptions. For each subject, the absolute difference between all possible pairs of intra-site and inter-site measures were computed. For instance, take the first row of the subject data.frame below, which has BWH as the “reference site”: the absolute differences between that row and the 6 non-BWH rows are computed for \(y\)——this process is repeated for all rows (and the results averaged) to arrive to the average inter-site difference for this subject; for the average intra-site differences the absolute difference of the 4 intrasite pairs are averaged.

df %>% 
  filter(subject == '01001') %>% 
  select(!where(is.numeric), ATROPOS_GM_AD) %>% 
  rename(y = ATROPOS_GM_AD)

These differences were averaged across all subjects resulting in the Mean Absolute Difference Between Sites (\(\overline{MAD_{B}}\)) and the Mean Absolute Difference Within Sites (\(\overline{MAD_{W}}\)). The test statistic is composed of their ratio:

\[ F = \frac{\overline{MAD_{B}}}{\overline{MAD_{W}}} \]

Key Findings: Permutation tests revealed significant site effects in that all 36 outcomes.

Data Inventory

Subjects 03001 and 03002 didn’t complete imaging at all 4 sites. Subject 03001 completed imaging at Penn and BWH only and 03002 completed imaging at BWH, Hopkins, and the NIH only.

Subject 02001 is missing missing DTI data from their NIH visit and is therefore excluded from this analysis.

Scans at site per subject

t_df <- df %>% 
   count(subject, site, .drop = FALSE) %>% 
  tidyr::pivot_wider(names_from = site, values_from = n) %>% 
  mutate(`Row Totals` = BWH + Hopkins + NIH + Penn) 
t_df <- t_df %>% 
  summarize(across(where(is.integer), sum)) %>% 
  bind_rows(t_df, .) %>% 
  mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_df

For Atropos (substracting failed segmentations)

Atropos segmentation failed for 6 images: 04001NIH01, 01003NIH01, 03002NIH02, 04003NIH01, 04001BWH02, and 03001BWH01.

atropos fail

#knitr::include_graphics(here::here(c('misc/atropos_fail.png', 'misc/atropos_success.png')))
#knitr::include_graphics("misc/atropos_fail.png")
#knitr::include_graphics("misc/atropos_success.png")
t_df <- df %>% 
  unite(sub, subject, site, scan, sep = '-') %>%
  filter(!sub %in% c('04001-NIH-01', '01003-NIH-01', '03002-NIH-02', '04003-NIH-01', '04001-BWH-02', '03001-BWH-01')) %>%
  separate(sub, c('subject', 'site', 'scan')) %>%
  mutate_if(is.character, as.factor) %>% 
  count(subject, site, .drop = FALSE) %>% 
  tidyr::pivot_wider(names_from = site, values_from = n) %>% 
  mutate(`Row Totals` = BWH + Hopkins + NIH + Penn)

t_df <- t_df %>% 
  summarize(across(where(is.integer), sum)) %>% 
  bind_rows(t_df, .) %>% 
  mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_df

Visualizations

Fractional Anisotropy

Densities

Densities are colored by site, while the black line represents the overall density aggregated across sites.

plot_densities <- function(var){
  df %>% 
    ggplot(aes_string(x=var, color='site')) + 
    geom_density() + 
    stat_density(aes_string(x = var), geom = "line", color = "black") +
    theme_bw() + 
    xlab(str_replace_all(var, "_", " "))
}

vars <-df %>% 
  select(ends_with('FA')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM FA

ATROPOS WM FA

FAST GM FA

FAST WM FA

FIRST THALAMUS FA

JLF GM FA

JLF WM FA

JLF THALAMUS FA

MIMOSA LESION FA

Boxplots

plot_avg_tensor_values <- function(tensor){
  df %>% 
    select(!is.numeric, ends_with(tensor)) %>% 
    pivot_longer(ends_with(tensor), names_to = 'seg', values_to = 'average') %>% 
    separate(seg, into = c('segmentation', 'roi', 'tensor')) %>% 
    unite(segmentation, segmentation, roi, sep = " ") %>% 
    ggplot(aes(x=site, y=average)) + 
    geom_boxplot(outlier.shape = NA) + 
    geom_jitter(aes(color = subject, shape=scan), alpha=0.6) +
    #facet_grid(site~.) + 
    facet_grid(segmentation ~ .) +
    coord_flip() +
    ggtitle(sprintf("Mean %s in ROI by site", tensor)) +
    theme_bw() + 
    # theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1)) +
    scale_x_discrete(expand = c(.00000001, .00000001)) +
    xlab("") + 
    ylab("") + 
    theme(strip.text.y = element_text(angle = 0))
}

plot_avg_tensor_values('FA')

Mean Diffusivity

Densities

vars <-df %>% 
  select(ends_with('MD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####", str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM MD

ATROPOS WM MD

FAST GM MD

FAST WM MD

FIRST THALAMUS MD

JLF GM MD

JLF WM MD

JLF THALAMUS MD

MIMOSA LESION MD

Boxplot

plot_avg_tensor_values('MD')

Radial Diffusivity

Densities

vars <-df %>% 
  select(ends_with('RD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM RD

ATROPOS WM RD

FAST GM RD

FAST WM RD

FIRST THALAMUS RD

JLF GM RD

JLF WM RD

JLF THALAMUS RD

MIMOSA LESION RD

Boxplot

plot_avg_tensor_values('RD')

Axial Diffusivity

Densities

vars <-df %>% 
  select(ends_with('AD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM AD

ATROPOS WM AD

FAST GM AD

FAST WM AD

FIRST THALAMUS AD

JLF GM AD

JLF WM AD

JLF THALAMUS AD

MIMOSA LESION AD

Boxplot

plot_avg_tensor_values('AD')

Testing of Site Effects

pseudo-F Ratio statistics: Intersite Variability / Intrasite Variability

ratio_stats = readRDS("results/avg_tensor_by_roi_wide_ratio_stat.rds")
null_dist = readRDS("results/avg_tensor_by_roi_wide_null.rds")
ratio_stats

Distribution of permuted pseudo-F stats:

null_dist %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=name, y=value)) + 
  geom_violin(draw_quantiles = c(0.95)) + 
  geom_point(aes(x=name, y=value), 
             data = ratio_stats %>% 
               pivot_longer(everything())) + 
  coord_flip()

P-value of pseudo-F stat:

test_ratio_stat = function(ratio.stats, null.dists){
  
  p.vals = c()
  for (col in colnames(ratio.stats)){
    ratio.stat = ratio.stats[, col]
    null.dist = null.dists[, col]
    
    percentile = ecdf(null.dist)
    p.vals = c(p.vals, percentile(ratio.stat))
  }
  names(p.vals) = colnames(ratio.stats)
  return(p.vals)
}

p = 1 - test_ratio_stat(ratio_stats, null_dist)
p_df <- data.frame(p_val = (10000*p + 1)/(1+10000)) 
p_df %>% 
  mutate(p_val_fdr = p.adjust(p_val, method = 'fdr'))

Mixed-Effect Models

models <- outcomes %>% 
  lapply(function(x) lmer(reformulate("(1|subject) + (1|site:subject)", x), data=df)) %>% 
  setNames(outcomes)
icc_df <- models %>% 
  lapply(performance::icc, by_group = TRUE) %>% 
  bind_rows(.id = "outcome") %>% 
  as.data.frame()
icc_df %>% 
  filter(Group == 'subject') %>% 
  ggplot(aes(x=outcome, y=ICC)) + 
  geom_point() + 
  coord_flip() + 
  ggtitle("ICC of Subject Random Effect") + 
  ylim(0,1)

icc_df %>% 
  filter(Group == 'site:subject') %>% 
  ggplot(aes(x=outcome, y=ICC)) + 
  geom_point() + 
  coord_flip() + 
  ggtitle("ICC of Site:Subject Random Effect") + 
  ylim(0,1)

Diagnostics

perf_mixed_model = function(model){
  
  list(summary = summary(model),
       perf = performance::icc(model, by_group = TRUE),
       shapiro.resid = shapiro.test(residuals(model)),
       shapiro.subject = shapiro.test(coef(model)$`subject`[,1]),
       `shapiro.site:subject` = shapiro.test(coef(model)$`site:subject`[,1]),
       
       qq.resid = car::qqPlot(residuals(model),dist="norm", main="Residuals"),
       qq.ranef.subject = car::qqPlot(ranef(model)$subject[,1],dist="norm", main="Rand. Eff. subject"),
       `qq.ranef.site:subject` = car::qqPlot(ranef(model)$`site:subject`[,1],dist="norm", main="Rand. Eff. site:subject"),
       res.v.fit = ggplot(mapping = aes(y = resid(model), x = fitted(model))) + geom_abline(intercept = 0, slope = 0) + geom_point() + geom_smooth() + theme_bw()
       )
}

FA

ATROPOS WM
perf_mixed_model(models$ATROPOS_WM_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -507.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2705 -0.4084 -0.0074  0.4260  1.8737 
## 
## Random effects:
##  Groups       Name        Variance   Std.Dev.
##  site:subject (Intercept) 0.00004437 0.006661
##  subject      (Intercept) 0.00042895 0.020711
##  Residual                 0.00001039 0.003223
## Number of obs: 74, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.387130   0.006349   60.97
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.092
## subject      | 0.887
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.97572, p-value = 0.1649
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.94623, p-value = 0.5962
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.9703, p-value = 0.3677
## 
## 
## $qq.resid
## 53 37 
## 51 36 
## 
## $qq.ranef.subject
## [1] 11  1
## 
## $`qq.ranef.site:subject`
## [1]  8 18
## 
## $res.v.fit

ATROPOS GM
perf_mixed_model(models$ATROPOS_GM_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -602.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1815 -0.4246  0.0087  0.2757  4.3686 
## 
## Random effects:
##  Groups       Name        Variance    Std.Dev.
##  site:subject (Intercept) 0.000005326 0.002308
##  subject      (Intercept) 0.000038020 0.006166
##  Residual                 0.000005805 0.002409
## Number of obs: 74, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.160945   0.001919   83.87
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.108
## subject      | 0.774
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.82302, p-value = 0.0000000553
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.98866, p-value = 0.9958
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.96679, p-value = 0.2835
## 
## 
## $qq.resid
## 51 52 
## 49 50 
## 
## $qq.ranef.subject
## [1] 9 2
## 
## $`qq.ranef.site:subject`
## [1]  8 18
## 
## $res.v.fit

FAST WM
perf_mixed_model(models$FAST_WM_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -557.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23611 -0.48150  0.01147  0.42822  2.17374 
## 
## Random effects:
##  Groups       Name        Variance   Std.Dev.
##  site:subject (Intercept) 0.00002158 0.004646
##  subject      (Intercept) 0.00057907 0.024064
##  Residual                 0.00001383 0.003719
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.403259   0.007307   55.19
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.035
## subject      | 0.942
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.97818, p-value = 0.1882
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.93861, p-value = 0.5043
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.97583, p-value = 0.5383
## 
## 
## $qq.resid
## [1] 32 31
## 
## $qq.ranef.subject
## [1]  1 11
## 
## $`qq.ranef.site:subject`
## [1] 24  8
## 
## $res.v.fit

FAST GM
perf_mixed_model(models$FAST_GM_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -639.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4185 -0.3794  0.0231  0.3556  4.4690 
## 
## Random effects:
##  Groups       Name        Variance    Std.Dev.
##  site:subject (Intercept) 0.000006098 0.002469
##  subject      (Intercept) 0.000059127 0.007689
##  Residual                 0.000006814 0.002610
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.173461   0.002371   73.15
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.085
## subject      | 0.821
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.84112, p-value = 0.00000008246
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.94055, p-value = 0.5269
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.98777, p-value = 0.9369
## 
## 
## $qq.resid
## [1] 51 52
## 
## $qq.ranef.subject
## [1] 9 2
## 
## $`qq.ranef.site:subject`
## [1]  8 18
## 
## $res.v.fit

JLF WM
perf_mixed_model(models$JLF_WM_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -570.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15403 -0.62253  0.05912  0.47991  1.65535 
## 
## Random effects:
##  Groups       Name        Variance    Std.Dev.
##  site:subject (Intercept) 0.000050221 0.007087
##  subject      (Intercept) 0.000468884 0.021654
##  Residual                 0.000006638 0.002576
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.426441   0.006635   64.27
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.096
## subject      | 0.892
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.99012, p-value = 0.8027
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.88581, p-value = 0.1232
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.95946, p-value = 0.1605
## 
## 
## $qq.resid
## [1] 53 54
## 
## $qq.ranef.subject
## [1] 11  2
## 
## $`qq.ranef.site:subject`
## [1]  8 24
## 
## $res.v.fit

JLF GM
perf_mixed_model(models$JLF_GM_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -671.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.32410 -0.42748 -0.05894  0.43779  2.37595 
## 
## Random effects:
##  Groups       Name        Variance    Std.Dev.
##  site:subject (Intercept) 0.000015087 0.003884
##  subject      (Intercept) 0.000041948 0.006477
##  Residual                 0.000002227 0.001492
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.149719   0.002058   72.75
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.255
## subject      | 0.708
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.96707, p-value = 0.03697
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.93611, p-value = 0.4759
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.96696, p-value = 0.2871
## 
## 
## $qq.resid
## [1] 78 77
## 
## $qq.ranef.subject
## [1] 9 2
## 
## $`qq.ranef.site:subject`
## [1] 28 27
## 
## $res.v.fit

FIRST THAL
perf_mixed_model(models$FIRST_THALAMUS_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -510.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.49646 -0.51505  0.04177  0.48047  1.81130 
## 
## Random effects:
##  Groups       Name        Variance   Std.Dev.
##  site:subject (Intercept) 0.00007717 0.008785
##  subject      (Intercept) 0.00036529 0.019113
##  Residual                 0.00002168 0.004656
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.309516   0.005958   51.95
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.166
## subject      | 0.787
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98419, p-value = 0.4289
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.9135, p-value = 0.2682
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.96969, p-value = 0.3518
## 
## 
## $qq.resid
## [1] 52 70
## 
## $qq.ranef.subject
## [1] 2 7
## 
## $`qq.ranef.site:subject`
## [1] 24 26
## 
## $res.v.fit

JLF THAL
perf_mixed_model(models$JLF_THALAMUS_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -464.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6797 -0.4081  0.0280  0.3874  3.4486 
## 
## Random effects:
##  Groups       Name        Variance   Std.Dev.
##  site:subject (Intercept) 0.00019846 0.014088
##  subject      (Intercept) 0.00023697 0.015394
##  Residual                 0.00003798 0.006163
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.319975   0.005212    61.4
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.419
## subject      | 0.501
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.90137, p-value = 0.00001392
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.90554, p-value = 0.2158
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.96625, p-value = 0.2722
## 
## 
## $qq.resid
## [1] 51 52
## 
## $qq.ranef.subject
## [1] 1 9
## 
## $`qq.ranef.site:subject`
## [1] 10  8
## 
## $res.v.fit

MIMOSA
perf_mixed_model(models$MIMOSA_LESION_FA)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -385.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7344 -0.3797  0.0021  0.3696  3.2023 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev.
##  site:subject (Intercept) 0.0003639 0.01908 
##  subject      (Intercept) 0.0013332 0.03651 
##  Residual                 0.0001144 0.01070 
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   0.3258     0.0115   28.34
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.201
## subject      | 0.736
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.93547, p-value = 0.0005552
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.97161, p-value = 0.9023
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.87229, p-value = 0.0003283
## 
## 
## $qq.resid
## [1] 24 23
## 
## $qq.ranef.subject
## [1] 4 5
## 
## $`qq.ranef.site:subject`
## [1] 14 33
## 
## $res.v.fit

MD

ATROPOS WM
perf_mixed_model(models$ATROPOS_WM_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1549.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5582 -0.3564  0.0353  0.4075  2.8736 
## 
## Random effects:
##  Groups       Name        Variance          Std.Dev.   
##  site:subject (Intercept) 0.000000000046271 0.000006802
##  subject      (Intercept) 0.000000000236783 0.000015388
##  Residual                 0.000000000004724 0.000002174
## Number of obs: 74, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000616523 0.000004775   129.1
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.161
## subject      | 0.823
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.81595, p-value = 0.00000003449
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.93972, p-value = 0.5172
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.95053, p-value = 0.0791
## 
## 
## $qq.resid
## 53 54 
## 51 52 
## 
## $qq.ranef.subject
## [1] 1 3
## 
## $`qq.ranef.site:subject`
## [1] 18 21
## 
## $res.v.fit

ATROPOS GM
perf_mixed_model(models$ATROPOS_GM_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1503.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.65330 -0.40529  0.02173  0.35056  2.27089 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000007871 0.000008872
##  subject      (Intercept) 0.00000000020654 0.000014371
##  Residual                 0.00000000001163 0.000003410
## Number of obs: 74, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000754518 0.000004583   164.6
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.265
## subject      | 0.696
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.95972, p-value = 0.0188
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.90487, p-value = 0.2118
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.97157, p-value = 0.403
## 
## 
## $qq.resid
## [1] 14 13
## 
## $qq.ranef.subject
## [1] 8 1
## 
## $`qq.ranef.site:subject`
## [1] 32 12
## 
## $res.v.fit

FAST WM
perf_mixed_model(models$FAST_WM_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1671.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5150 -0.4235 -0.0098  0.4002  2.8842 
## 
## Random effects:
##  Groups       Name        Variance          Std.Dev.   
##  site:subject (Intercept) 0.000000000063511 0.000007969
##  subject      (Intercept) 0.000000000262843 0.000016212
##  Residual                 0.000000000005098 0.000002258
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000611772 0.000005061   120.9
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.192
## subject      | 0.793
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.86709, p-value = 0.000000636
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.93033, p-value = 0.4142
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.91937, p-value = 0.007356
## 
## 
## $qq.resid
## [1] 53 54
## 
## $qq.ranef.subject
## [1] 1 3
## 
## $`qq.ranef.site:subject`
## [1] 18 21
## 
## $res.v.fit

FAST GM
perf_mixed_model(models$FAST_GM_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1644.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.57814 -0.48669 -0.01082  0.39656  2.24634 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000006816 0.000008256
##  subject      (Intercept) 0.00000000016182 0.000012721
##  Residual                 0.00000000001044 0.000003231
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000736558 0.000004076   180.7
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.284
## subject      | 0.673
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.96328, p-value = 0.02141
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.91875, p-value = 0.3084
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.98161, p-value = 0.7491
## 
## 
## $qq.resid
## [1] 14 13
## 
## $qq.ranef.subject
## [1] 1 8
## 
## $`qq.ranef.site:subject`
## [1] 18 22
## 
## $res.v.fit

JLF WM
perf_mixed_model(models$JLF_WM_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1655
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.10744 -0.47068  0.05638  0.40899  2.08380 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000002550 0.000005050
##  subject      (Intercept) 0.00000000028330 0.000016832
##  Residual                 0.00000000001291 0.000003593
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000649895 0.000005156     126
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.079
## subject      | 0.881
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98922, p-value = 0.7456
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.91596, p-value = 0.2864
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.97438, p-value = 0.4894
## 
## 
## $qq.resid
## [1] 53 35
## 
## $qq.ranef.subject
## [1] 1 3
## 
## $`qq.ranef.site:subject`
## [1] 18 22
## 
## $res.v.fit

JLF GM
perf_mixed_model(models$JLF_GM_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1624.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6447 -0.3965  0.0703  0.4914  1.8407 
## 
## Random effects:
##  Groups       Name        Variance          Std.Dev.   
##  site:subject (Intercept) 0.000000000121841 0.000011038
##  subject      (Intercept) 0.000000000461916 0.000021492
##  Residual                 0.000000000008956 0.000002993
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000802476 0.000006729   119.3
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.206
## subject      | 0.779
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.99124, p-value = 0.8673
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.92387, p-value = 0.3522
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.97873, p-value = 0.6423
## 
## 
## $qq.resid
## [1] 58 57
## 
## $qq.ranef.subject
## [1] 8 1
## 
## $`qq.ranef.site:subject`
## [1] 18  6
## 
## $res.v.fit

FIRST THAL
perf_mixed_model(models$FIRST_THALAMUS_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1585.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.21713 -0.55661  0.09266  0.56833  2.12975 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000006565 0.000008103
##  subject      (Intercept) 0.00000000021540 0.000014676
##  Residual                 0.00000000003787 0.000006154
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000664029 0.000004667   142.3
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.206
## subject      | 0.675
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98821, p-value = 0.6786
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.92793, p-value = 0.3903
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.96794, p-value = 0.309
## 
## 
## $qq.resid
## [1] 21 35
## 
## $qq.ranef.subject
## [1] 3 7
## 
## $`qq.ranef.site:subject`
## [1] 22 24
## 
## $res.v.fit

JLF THAL
perf_mixed_model(models$JLF_THALAMUS_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1574.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.07380 -0.40130  0.02501  0.52883  1.57687 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000009030 0.000009503
##  subject      (Intercept) 0.00000000030257 0.000017394
##  Residual                 0.00000000003792 0.000006158
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000651264 0.000005509   118.2
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.210
## subject      | 0.702
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.97973, p-value = 0.2348
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.93132, p-value = 0.4243
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.98705, p-value = 0.9206
## 
## 
## $qq.resid
## [1] 51 21
## 
## $qq.ranef.subject
## [1] 1 3
## 
## $`qq.ranef.site:subject`
## [1] 22 29
## 
## $res.v.fit

MIMOSA
perf_mixed_model(models$MIMOSA_LESION_MD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_MD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1352.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.90417 -0.44509 -0.02298  0.38772  1.58995 
## 
## Random effects:
##  Groups       Name        Variance       Std.Dev.  
##  site:subject (Intercept) 0.000000002872 0.00005359
##  subject      (Intercept) 0.000000003487 0.00005905
##  Residual                 0.000000000462 0.00002149
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00091264 0.00001993   45.79
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.421
## subject      | 0.511
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.99023, p-value = 0.8093
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.94465, p-value = 0.5766
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.9748, p-value = 0.5034
## 
## 
## $qq.resid
## [1] 24 19
## 
## $qq.ranef.subject
## [1] 6 7
## 
## $`qq.ranef.site:subject`
## [1] 33 14
## 
## $res.v.fit

RD

ATROPOS WM
perf_mixed_model(models$ATROPOS_WM_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1569.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.73497 -0.46998 -0.01091  0.49515  1.76890 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000002355 0.000004853
##  subject      (Intercept) 0.00000000040997 0.000020248
##  Residual                 0.00000000000382 0.000001954
## Number of obs: 74, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00047857 0.00000616   77.69
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.054
## subject      | 0.937
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.9918, p-value = 0.9183
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.91991, p-value = 0.3179
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.93204, p-value = 0.01878
## 
## 
## $qq.resid
## 38 37 
## 37 36 
## 
## $qq.ranef.subject
## [1]  1 11
## 
## $`qq.ranef.site:subject`
## [1] 21 30
## 
## $res.v.fit

ATROPOS GM
perf_mixed_model(models$ATROPOS_GM_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1502.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.64340 -0.40411  0.02559  0.35378  2.07830 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000006018 0.000007757
##  subject      (Intercept) 0.00000000024306 0.000015590
##  Residual                 0.00000000001381 0.000003716
## Number of obs: 74, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000693166 0.000004886   141.9
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.190
## subject      | 0.767
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.92446, p-value = 0.0002792
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.9163, p-value = 0.289
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.98143, p-value = 0.7422
## 
## 
## $qq.resid
## 14 51 
## 14 49 
## 
## $qq.ranef.subject
## [1] 8 1
## 
## $`qq.ranef.site:subject`
## [1] 32 12
## 
## $res.v.fit

FAST WM
perf_mixed_model(models$FAST_WM_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1675.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78643 -0.43514  0.02686  0.39648  1.86959 
## 
## Random effects:
##  Groups       Name        Variance          Std.Dev.   
##  site:subject (Intercept) 0.000000000038492 0.000006204
##  subject      (Intercept) 0.000000000492182 0.000022185
##  Residual                 0.000000000005604 0.000002367
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000468050 0.000006769   69.14
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.072
## subject      | 0.918
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.99336, p-value = 0.9574
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.91766, p-value = 0.2996
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.95814, p-value = 0.1447
## 
## 
## $qq.resid
## [1] 38 37
## 
## $qq.ranef.subject
## [1]  1 11
## 
## $`qq.ranef.site:subject`
## [1] 21 18
## 
## $res.v.fit

FAST GM
perf_mixed_model(models$FAST_GM_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1642.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6925 -0.3926  0.0002  0.3444  2.4629 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000005085 0.000007131
##  subject      (Intercept) 0.00000000021045 0.000014507
##  Residual                 0.00000000001245 0.000003528
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000671595 0.000004541   147.9
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.186
## subject      | 0.769
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.90592, p-value = 0.00002186
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.95301, p-value = 0.6826
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.97631, p-value = 0.555
## 
## 
## $qq.resid
## [1] 51 14
## 
## $qq.ranef.subject
## [1] 1 8
## 
## $`qq.ranef.site:subject`
## [1] 32 18
## 
## $res.v.fit

JLF WM
perf_mixed_model(models$JLF_WM_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1682.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12628 -0.51321 -0.01984  0.35191  2.97520 
## 
## Random effects:
##  Groups       Name        Variance          Std.Dev.   
##  site:subject (Intercept) 0.000000000019247 0.000004387
##  subject      (Intercept) 0.000000000469674 0.000021672
##  Residual                 0.000000000007354 0.000002712
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00048768 0.00000658   74.12
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.039
## subject      | 0.946
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.9634, p-value = 0.02177
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.9182, p-value = 0.3039
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.91105, p-value = 0.004072
## 
## 
## $qq.resid
## [1] 35 36
## 
## $qq.ranef.subject
## [1]  1 11
## 
## $`qq.ranef.site:subject`
## [1] 22 28
## 
## $res.v.fit

JLF GM
perf_mixed_model(models$JLF_GM_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1636.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.54869 -0.36475  0.00208  0.50417  1.80211 
## 
## Random effects:
##  Groups       Name        Variance          Std.Dev.   
##  site:subject (Intercept) 0.000000000093489 0.000009669
##  subject      (Intercept) 0.000000000517566 0.000022750
##  Residual                 0.000000000007899 0.000002811
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000744510 0.000007042   105.7
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.151
## subject      | 0.836
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98258, p-value = 0.3476
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.93028, p-value = 0.4137
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.98775, p-value = 0.9365
## 
## 
## $qq.resid
## [1] 58 57
## 
## $qq.ranef.subject
## [1] 8 1
## 
## $`qq.ranef.site:subject`
## [1]  6 18
## 
## $res.v.fit

FIRST THAL
perf_mixed_model(models$FIRST_THALAMUS_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1600.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64679 -0.48633 -0.05519  0.53507  2.06990 
## 
## Random effects:
##  Groups       Name        Variance         Std.Dev.   
##  site:subject (Intercept) 0.00000000005445 0.000007379
##  subject      (Intercept) 0.00000000022682 0.000015060
##  Residual                 0.00000000002980 0.000005459
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000554717 0.000004736   117.1
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.175
## subject      | 0.729
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98806, p-value = 0.6687
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.98161, p-value = 0.9745
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.93314, p-value = 0.02042
## 
## 
## $qq.resid
## [1] 35 21
## 
## $qq.ranef.subject
## [1] 7 2
## 
## $`qq.ranef.site:subject`
## [1] 22 29
## 
## $res.v.fit

JLF THAL
perf_mixed_model(models$JLF_THALAMUS_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1564.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1758 -0.3921  0.0092  0.4768  2.4299 
## 
## Random effects:
##  Groups       Name        Variance        Std.Dev.   
##  site:subject (Intercept) 0.0000000001390 0.000011792
##  subject      (Intercept) 0.0000000002651 0.000016283
##  Residual                 0.0000000000382 0.000006181
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000539595 0.000005311   101.6
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.314
## subject      | 0.599
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.94219, p-value = 0.001264
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.81838, p-value = 0.01642
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.97448, p-value = 0.4927
## 
## 
## $qq.resid
## [1] 51 52
## 
## $qq.ranef.subject
## [1] 1 6
## 
## $`qq.ranef.site:subject`
## [1] 22 10
## 
## $res.v.fit

MIMOSA
perf_mixed_model(models$MIMOSA_LESION_RD)

## $summary
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_RD ~ (1 | subject) + (1 | site:subject)
##    Data: df
## 
## REML criterion at convergence: -1352.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.36040 -0.45821 -0.06635  0.43847  1.65410 
## 
## Random effects:
##  Groups       Name        Variance        Std.Dev.  
##  site:subject (Intercept) 0.0000000028530 0.00005341
##  subject      (Intercept) 0.0000000031759 0.00005635
##  Residual                 0.0000000004712 0.00002171
## Number of obs: 80, groups:  site:subject, 40; subject, 11
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 0.0007483  0.0000192   38.98
## 
## $perf
## # ICC by Group
## 
## Group        |   ICC
## --------------------
## site:subject | 0.439
## subject      | 0.489
## 
## $shapiro.resid
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.97547, p-value = 0.1267
## 
## 
## $shapiro.subject
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$subject[, 1]
## W = 0.9102, p-value = 0.2452
## 
## 
## $`shapiro.site:subject`
## 
##  Shapiro-Wilk normality test
## 
## data:  coef(model)$`site:subject`[, 1]
## W = 0.99057, p-value = 0.9811
## 
## 
## $qq.resid
## [1] 24 19
## 
## $qq.ranef.subject
## [1] 11  4
## 
## $`qq.ranef.site:subject`
## [1] 14 33
## 
## $res.v.fit

Repeated Measures ANOVA

library(rstatix)
find_icc <- function(data){
  
  irr::icc(data,
    model = "twoway", 
    type = "agreement",
    unit = "single")
  
}

run_anova <- function (var){
  df %>% 
    filter(scan == '01') %>% 
    anova_test(dv = var, wid = subject, within = site) %>% 
    get_anova_table()
}

anovas <- outcomes %>% 
  lapply(run_anova) %>% 
  setNames(outcomes)


data.frame(measure = names(anovas), 
           F_stat = anovas %>% 
             purrr::map(~ .x$F) %>% 
             purrr::reduce(c),
           p = anovas %>% 
             purrr::map(~ .x$p) %>% 
             purrr::reduce(c))
pair_df_list <- df %>% 
  pivot_longer(where(is.numeric)) %>% 
  split(list(.$site, .$name))

# compute names for each pair_df
df_names <- pair_df_list %>% 
  purrr::map(~ .x %>% 
               transmute(df_name = paste(site, name, sep="-")) %>% 
               pull() %>% 
               unique()
             ) %>% 
  purrr::reduce(c)
# finish transformation
pair_df_list <- pair_df_list %>% 
  purrr::map(~ .x %>% 
               pivot_wider(names_from = c('name', 'scan')) %>% 
               select(where(is.numeric))
               ) %>% 
  setNames(df_names)

icc_df <- pair_df_list %>% 
  purrr::map(find_icc) %>% 
  purrr::imap(~ data.frame(outcome = .y, 
                          icc = .x$value, 
                          lwr = .x$lbound,
                          upr = .x$ubound,
                          p_val = .x$p.value)) %>% 
  dplyr::bind_rows(.id = 'outcome') %>% 
  separate(outcome, into = c('site', 'outcome'), sep='-')
icc_df